Spatial characterization of the effect of age and sex on macular layer thicknesses and foveal pit morphology

Characterizing the effect of age and sex on macular retinal layer thicknesses and foveal pit morphology is crucial to differentiating between natural and disease-related changes. We applied advanced image analysis techniques to optical coherence tomography (OCT) to: 1) enhance the spatial description of age and sex effects, and 2) create a detailed open database of normative retinal layer thickness maps and foveal pit shapes. The maculae of 444 healthy subjects (age range 21–88) were imaged with OCT. Using computational spatial data analysis, thickness maps were obtained for retinal layers and averaged into 400 (20 x 20) sectors. Additionally, the geometry of the foveal pit was radially analyzed by computing the central foveal thickness, rim height, rim radius, and mean slope. The effect of age and sex on these parameters was analyzed with multiple regression mixed-effects models. We observed that the overall age-related decrease of the total retinal thickness (TRT) (-1.1% per 10 years) was mainly driven by the ganglion cell-inner plexiform layer (GCIPL) (-2.4% per 10 years). Both TRT and GCIPL thinning patterns were homogeneous across the macula when using percentual measurements. Although the male retina was 4.1 μm thicker on average, the greatest differences were mainly present for the inner retinal layers in the inner macular ring (up to 4% higher TRT than in the central macula). There was an age-related decrease in the rim height (1.0% per 10 years) and males had a higher rim height, shorter rim radius, and steeper mean slope. Importantly, the radial analysis revealed that these changes are present and relatively uniform across angular directions. These findings demonstrate the capacity of advanced analysis of OCT images to enhance the description of the macula. This, together with the created dataset, could aid the development of more accurate diagnosis models for macular pathologies.


Introduction
The macula is an approximately 5.5 mm diameter oval pigmented structure of the retina located near its center that is responsible for central high-resolution vision and color perception [1]. Compared to the peripheral retina, the cellular arrangement of the macular retina is unique: it has a high density of photoreceptors and a characteristic centripetal displacement of inner cellular layers that form a concave depression in its center called the fovea, upon which the visual axis is fixed. The macula has a crucial significance in human vision, indeed, pathologies affecting its integrity (e.g., age-related macular degeneration) can lead to severe visual impairment [2].
Optical coherence tomography (OCT)-a non-invasive imaging technique capable of obtaining micrometer resolution images of the retinal structure-has greatly facilitated the investigation of disease-related changes in the structural features of the macula [3]. From OCT images, computational image processing techniques have revealed alterations in macular layer thicknesses [4][5][6] as well as foveal pit morphology [7,8] due to ophthalmological and neurological conditions. Importantly, studies on healthy cohorts have shown that demographical factors such as age, sex, and ethnicity influence the structural parameters of the macula obtained from OCT images. For instance, age-related thinning of the inner retinal layers has been reported in a review [9]. Sex-differences in both thickness and foveal pit geometry have also been detected [10][11][12]. These findings evidence the importance of accurately characterizing the variation of the retinal structure in healthy populations so that robust conclusions can be reached from clinical studies. In this regard, current knowledge on the matter is limited by several aspects. First, most studies analyzing macular thickness have relied on the standard Early Treatment Diabetic Retinopathy Study (ETDRS) sectorization ( Fig 1B) [13][14][15][16][17][18][19][20][21][22]. This limits the description of the macular region to the average measurements of nine sectors within a 6-mm circle centered on the fovea. This choice might undermine the ability to capture more detailed spatial variation patterns and recent studies have begun to use smaller than ETDRS sectors. For instance, in [23] authors used 64 sectors (8 x 8 square grid) and found specific clustered spatial patterns of age changes. Similarly, the same 8 x 8 grid was used to establish a normative database of macular thickness in [24]. Additionally, a radial grid with 61 sectors has also been used to examine age changes in retinal thickness [25].
On the other hand, the limited published data available would suggest that the effect of both sex and age on the foveal pit geometry has only been partially studied to date. In the case of sex, after early work on the topic [10,26], the work of Scheibe et al. was the first relatively large study reporting clear sex differences in foveal pit morphology [11]. More recently, sex differences in foveal curvature were also found in a large study using the UK-Biobank dataset [27]. As for age, this factor has been less explored and large studies have focused mainly on the foveal slope [27][28][29]. Importantly, only a few works have investigated the foveal shape across multiple angular directions [11,30].
In light of this gap, we studied both retinal thickness and foveal pit morphology in a large sample of 444 healthy subjects. Following a multiscale approach, we evaluated age-related changes and sex differences on retinal thickness using both the regular ETDRS and a 20 x 20 square grid sectorizations. With the latter, we aimed to provide a more detailed normative database by using smaller sectors than those existing in the literature. Similarly, we assessed the effect of such demographic factors on the central thickness, rim height, radius, and slope of the foveal pit measured for the whole macula as well as for 24 angular directions individually. With this characterization we extend existing knowledge by examining changes in various aspects of the fovea with high spatial detail.

Participants
A total of 444 healthy subjects (855 eyes) were included in the study. Subjects were recruited at the Ophthalmology and Neurology Departments of the Cruces University Hospital (Barakaldo, Spain). The mean age of participants was 54.9 ± 12.7 years (ranging from 21 to 88) and 63% were females (see Table 1 for a summary and S1 Fig for the detailed age distribution). All subjects were Caucasian. Before inclusion in the study, all participants underwent a screening process that consisted of an ophthalmological examination and a comprehensive questionnaire on neurological, systemic, and eye-related diseases. We excluded subjects with a history  of heavy smoking (>20 cigarettes/day), heavy alcohol use (>4 drinks/day for men or >3 drinks/day for women), diagnosis of any type or grade of diabetes, uncontrolled or resistant elevated blood pressure, obesity (body mass index > 30), history of consumption of drugs or medications known to induce retinal toxicity, chronic inflammatory systemic diseases, history of traumatic brain injury, or neurological diseases. Additionally, we excluded subjects with spherical equivalent refractive error > 4.00 diopters or < -4.00 diopters, >3.00 diopters of astigmatism, or any other ocular condition potentially affecting OCT measures, as detailed in the OSCAR-IB consensus criteria [31]. In cases where only one of the eyes of a participant was excluded, the other eye was included. Following the tenets of the Declaration of Helsinki, all participants gave written informed consent prior to their participation. The study was approved by the Institutional Ethics Committee of OSI Ezkerraldea-Enkarterri-Cruces (Barakaldo, Spain).

Image acquisition and processing
The eyes of all subjects were imaged using a Spectralis spectral domain OCT scanner (Heidelberg Engineering, Heidelberg, Germany) following a macular raster acquisition protocol with a 20º field of view, 25 horizontal B-scans, and 512 A-scans per B-scan. For each final Bscan, a total of 49 B-scans were averaged. No pupil dilation was employed and default keratometry values were used. Additionally, both eyes of a subset of 12 subjects were imaged for a second time, covering the same macular area but following a higher resolution protocol with 97 B-Scans and 1024 A-scans. All other parameters of the second scans were the same. This dataset was used to investigate the impact of different acquisition protocols in a sensitivity analysis. Ocular magnification was corrected by the built-in Spectralis software [32]. All images were segmented with Heidelberg Eye Explorer 1.9.10.0 software. Segmentations were reviewed by three specialists (A.M., S.T. and I.G.) and evident errors within the 3 mm radius macular region were manually corrected. The images were loaded into MATLAB 2020b (The Mathworks Inc., Natick, MA, United States) and subsequently analyzed using the open-source RETIMAT Toolbox (https://github.com/drombas/retimat). All scans were aligned by automatically locating the foveal center as the minimum value of a smoothed total retinal thickness (TRT) map [33]. Left eyes were flipped to match right eyes. From the segmentation data, two macular features were studied: retinal layer thicknesses and foveal pit morphology. The retinal layers depicted in Fig 1A were studied: TRT, retinal nerve fiber layer (RNFL), ganglion cell-inner plexiform layer (GCIPL), inner nuclear layer (INL), outer nuclear and plexiform layer (ONPL), and external limiting membrane-Bruch's membrane complex (ELM-BM). Point thickness values were interpolated to a 300 x 300 regular grid, with a spacing of 0.1 mm. Then, these values were averaged following three fovea-centered sectorizations: whole macular region (i.e., the 3 mm radius circular region), ETDRS (Fig 1B), and a 6 x 6 mm sectorization with 20 x 20 square sectors ( Fig 1C). In this latter sectorization, regions outside of the 3 mm fovea-centered circle were excluded. Similarly, the thicknesses of the RNFL, GCIPL, and INL layers were not analyzed in the centermost sectors (i.e., the central 1.2 x 1.2 mm region). This is because the thicknesses of the inner layers in the central foveal region are almost zero and any results including these sectors could be highly biased by segmentation errors.
Regarding foveal pit morphology, TRT point values were interpolated to a radial grid with 24 angular directions, 2 mm radius, and 0.1 mm spacing. First, the central foveal thickness (CFT) was estimated as the TRT at the foveal center. LOESS smoothing (span = 50%) was then separately applied to the TRT profile for each radial direction to reduce the ripple before computing the rim height as the point of maximum TRT. Based on the rim height, two additional parameters were derived: rim radius i.e., the lateral distance from the foveal center to the rim, and mean slope i.e., the mean first derivative of the unsmoothed TRT between the foveal center and the rim (Fig 1D). Except for the CFT, which does not vary radially, all parameters were computed for both the whole macula (i.e., averaging all 24 directions) and for every single direction. S2 Fig illustrates this process.

Data analysis
The effects of sex and age on thickness maps and foveal pit parameters were studied in both absolute and percentual units by means of multivariate regression. The models included fixed terms for age and sex (with females as the reference). To adjust for differences in ocular shape, a fixed term for the scan focus variable was also included. This variable is estimated by the Spectralis scanner while focusing the image and accounts for the refraction error of each eye, which is known to influence retinal measurements. Importantly, a mixed-effects model variant with a random intercept (γ subject ) was used to account for the inter-eye correlation. As a first model selection step, two models with linear (Eq 1) and quadratic (Eq 2) age effects were fitted and the one with the minimum Akaike information criteria was chosen.
To equivalently report results from linear and quadratic models, a single combined age effect coefficient was estimated as the mean yearly change between the age of 40 and 80. For each coefficient, a 95% confidence interval (CI) and a p-value were calculated. When the selected model was quadratic, a single p-value was computed for the combined linear + quadratic effect of age by using an F-test to compare a reference model without any age term with the quadratic model. Age coefficients were transformed into % values by dividing the estimates by the average parameter value in the youngest group (age < = 40, including n = 51 from the 444 subjects). Likewise, sex coefficients were divided by the estimate for females in the youngest group (age < = 40, n = 37) to obtain % values. The significance level was set to 0.05. Results for the whole macular region were adjusted for multiple comparisons based on the Holm-Bonferroni correction [34]. For both ETDRS and the high-resolution sectorizations, a correction based on the false discovery rate [35] was employed due to the high number of sectors and the potential statistical dependence between tests. The overall fit of the model was assessed by the marginal R-squared, which measures the variance explained solely by the fixed terms. The age coefficient was reported as changes per 10 years as annual changes were found to be very small.
Using the subset of subjects imaged twice, we compared the measurements obtained by the standard protocol (25 B-scans, 512 A-scans) with the high-resolution protocol (97 B-scans, 1024 A-scans). For every parameter, a mixed-effects model linear regression was fitted with a fixed term for the bias introduced by using the standard protocol (β bias ) and a random intercept (γ subject ) to account for inter-eye correlation: A p-value as well as a 95% confidence interval (CI) were computed for the bias term. Finally, the estimated bias was divided by the intercept to obtain a percentual bias.

Effects of age and sex on macular thickness
Fig 2 shows the mean thickness as a function of age on each layer for the whole macular region. Table 2 reports the regression coefficient estimates along with the corresponding p-values and the R 2 of the model. All layers were found to become significantly thinner with age except for the RNFL (with a non-significant thickening and a high dispersion) and the ELM-BM. The thinning effect was stronger for the TRT, GCIPL, and INL. From these, the GCIPL presented the highest percentual loss (-2.4 [-2.9, -1.9] % every 10 years) and was responsible for most of the reduction in the TRT (-1.1 [-1.3, -0.8] % every 10 years).
On the other hand, male retinas were found to be 4.1 [1.8, 6.5] μm thicker on average ( Table 2). This was a combined effect of differences in all individual layers, although the INL, ONPL, and ELM-BM were the only layers that showed statistically significant differences between males and females. When the spatial distribution of thickness changes in relation to age and sex on the 20 x 20 grid was analyzed (Fig 3), we observed that the age-related thinning of the TRT, GCIPL and ONPL was relatively homogeneous across the macular area. For the INL, however, the agerelated reduction in thickness was more localized in the outer ring of the macula (perifoveal region). In the case of the RNFL maps, there was a mild thickening which was especially prominent in the temporal sector. Finally, the changes in the ELM-BM were minor and not significant except for the foveal central region where a slight-moderate thinning was observed. The regions in which a quadratic model was selected are depicted in S3 Fig. With regards to sex, differences in the TRT were prominent in the central region (radius<1.5 mm), with males having up to a 4% thicker retina. As for the individual layers, the sex differences in the RNFL, GCIPL, and INL layers were more obvious in the inner ring (0.5 mm � radius � 1.5 mm), where thickness values in males were also significantly higher than in females. At larger radii, these differences diminished and even reversed. Sex differences in the ONPL and ELM-BM layers were less pronounced and more homogeneous in percentage terms. However, they remained significant (males > females), predominantly for the outer ring (perifovea) in the ELM-BM and for the inner ring (parafovea) in the ONPL. For completeness, full regression coefficients including linear, quadratic and scan focus terms as well as correspondent results for ETDRS sectorization are reported as S1 Appendix.

Effects of age and sex on foveal pit morphology
Both sex and age influenced the foveal pit morphology (Table 3 and  Percentage differences were greatest for the mean slope (6.2 [2.8, 9.6] %).
The observed age and sex differences were also found when the foveal morphology was studied radially (Fig 4). More specifically, the estimated reduction in rim height was present in all directions and ranged from 0.6% to 0.8% per 10 years. The rim radius was found to decrease in all directions, although the effect was stronger in the inferior and nasal sectors, with a 1%

PLOS ONE
The effect of age and sex on the macula decrease per 10 years. Conversely, sex-related differences were clearly present in every angular direction. The differences were relatively uniformly distributed except for the rim radius, which showed slightly larger differences for the superior and inferior sectors.

Sensitivity analysis
The estimated bias introduced in each average thickness and foveal parameter due to the resolution of the acquisition protocol is shown in Table 4. Corresponding results for ETDRS and 20 x 20 sectorizations can be found in S1 Appendix and S5 Fig, respectively. The results highlight that the bias is small for average macular thickness (< 1%) but is not negligible for central thickness (3.02%) and foveal slope (-6.57%). Importantly, when using smaller sectors this bias is below 5% overall with the centermost and outer regions being the most affected.

PLOS ONE
The effect of age and sex on the macula

Discussion
In this study, we evaluated the influence of age and sex on the structure of the retina by testing their relationship with finely sectorized thickness maps of macular layers and foveal pit morphology metrics extracted from OCT images of 444 healthy subjects (855 eyes). We first observed that in relation to age there is a homogeneous TRT thinning of 1.1% per 10 years, driven mainly by the GCIPL (-2.4%) followed by the INL (-1.3%) and ONPL (-0.7%). The rim height decreased significantly with age at a rate of 1.0% every 10 years. Second, we found that on average male retinas were 4.1 μm thicker, with more prominent differences in the central region (radius<1.5 mm). Males had a larger central CFT, higher rim height, shorter rim radius and a steeper mean slope. The observed age-related TRT decline portrays the retina as an evolving structure, a finding that is well supported in the literature [9,15,[18][19][20][21]. In particular, one review of the literature described the age-related thinning pattern as spatially-dependent, with an unchanged or thickened central retina, and a maximum thickness loss in the parafoveal region [9]. In the present study, we focused on the percentage loss of TRT and found a relatively uniform thinning effect for eccentricities larger than 0.5 mm. This suggests that, except for the central macula, differences in the absolute thinning rate can be explained by differences in baseline TRT thickness.
Of greater interest than the TRT analysis, however, is the determination of the individual layers (with specific cellular architectures and spatial distributions) that drive the thinning effect. In this regard, histological studies have observed an age-related decrease in the number of fibers that comprise the optic nerve and the RNFL [36,37]. In contrast, we found that the RNFL remained unchanged or even became thicker in some temporal regions. This discrepancy is also present in the literature, with different studies reporting thinning [22,38], no effect [18,19], or thickening [20,39]. These mixed results might be explained by the thinness of the RNFL in relation to the axial resolution, which makes thickness measurements more prone to segmentation errors. Therefore, using peripapillary instead of macular OCT seems more appropriate to assess the RNFL.
As for the GCIPL, the thinning we observed is supported by histological evidence that points to an age-related decrease in ganglion cell density [37]. Similarly, the majority of previous studies reported a thinning of both the ganglion cell layer (GCL) [22,38,40] and the inner plexiform layer (IPL) [17,20,21,25]. We found the GCIPL to be the layer most affected

PLOS ONE
The effect of age and sex on the macula by age, which may indicate that it is particularly sensitive to aging. Additionally, we measured a consistent GCIPL percentual loss across the macula, except for the central region. These results are in line with previous work [41,42] and suggest that regional differences in absolute thickness loss-such as the accentuated parafoveal thinning [25]-reflect differences in absolute GCL thickness and not a spatially-dependent predisposition to an age-related decline. As for the INL, a thinning effect has been previously reported [18,20,21,25]. In a similar vein, our results describe the INL as the layer with the second most important age decline and a prevalent thinning of the outer regions. Finally, and in contrast to the thickening of the outer retinal layers described previously [9], we measured a thinning of both the ONPL and the ELM-BM, a finding more in line with recent studies [21,25,39]. Interestingly, we found two different thinning patterns: a prevalent and uniform ONPL thinning of outer regions, and a highly localized central region ELM-BM thinning. The different cellular configurations of these layers might explain these characteristic patterns.
Although the repeatability of Spectralis has been reported to be very good [43], it is important to note that the observed yearly changes are, for the most part, small when compared to the coefficient of variation. For instance, the coefficient of variation of Spectralis TRT measurements (with eye-tracking mode) is up to 0.86% [44], which would correspond to a decade change solely due to age. This underscores two points: 1) natural age changes would have a relatively small impact on longitudinal studies with a regular follow up (e.g., < 5 years) and 2), despite its high repeatability, OCT requires large sample sizes and groupwise statistical analysis to detect small changes in the retina.
Regarding sex differences, a thicker retina in males has been previously reported in both adults [13][14][15]18] and children [45]. In all studies, differences were higher for the inner macular ring and diminished for the outer ring, a pattern also observed in our analysis. In addition, our analysis revealed spatially localized differences in all layers, which suggests that sex plays a role in the entire retina. As further evidence of this, layer-specific differences have been described previously. For instance, some studies measured a thicker RNFL in males [17,18,21,24], while others observed it to be thicker in females [19,22,38]. In addition, a greater thickness in males has been found for the GCL [16,18,21,24], inner plexiform layer [17,18,21] and INL [19]. As for the outer retinal layers, sex differences have also been reported [17,19].
A possible explanation for this thicker retina in males could be a systematic macroscopic difference. In fact, a clear finding from magnetic resonance imaging studies is that males have a larger ocular globe [46] and, therefore, one could hypothesize that larger eyes have a thicker retina. Contrary to this, a negative relationship between axial length and retinal thickness has been observed [47], although this might be confounded by the ocular effect of ocular biometry on lateral image scaling [48,49] (i.e., the same field of view corresponds to a larger region in larger eyes).
Similarly, a sex bias in lateral image scale estimation could explain why sex differences in inner layer thicknesses are localized in the inner ring. Although lateral scaling is usually adjusted for axial length differences [10,21], scanners typically assume a nominal corneal curvature value for both males and females (e.g. Spectralis uses a 7.7 mm [32]). This assumption might lead to an overestimation of the lateral scaling in females, as corneal curvature is both positively related with lateral scaling and smaller in females [50]. This, in turn, could displace the maxima of the TRT, GCIPL and INL thickness profiles in females to larger eccentricities, thereby resulting in the observed pattern.
With respect to the effect of age on the foveal pit, neither previous studies [26,30] nor our data showed a clear effect on the CFT (S1 Table). This is likely to be related to the absence of inner retinal layers in the central region, which are more prone to becoming thinner [9].
Regarding the rim height, although smaller previous studies did not find a statistically significant effect [26,30] we observed a clear thinning of the rim, which is in line with the known TRT loss. In addition, we also observed a decrease in the rim radius, which might be a consequence of a flattening of the rim due to the rim height decrease. This, however, did not remain statistically significant after correction. The high inter-subject variability of the mean slope resulted in high uncertainty in the estimates which is likely responsible for the mixed results in the literature [12,28,29].
Although existing studies describing sex differences in the foveal pit used different parameter definitions and mathematical models, one finding is consistent between our results and most of the previous work: a broader and shallower pit in females [10][11][12]27, 51] (S2 Table). As with thickness differences, the lateral scale estimation bias introduced by ocular magnification might lead to an overestimation of the lateral scaling in females and the observed differences in slope and radius. As for height measurements, our study confirms previous findings of greater CFT and rim height in males [11]. Considering that these foveal pit metrics are effectively thickness measurements, it is likely that these differences are simply a consequence of a higher TRT in males.
In addition to the overall description, the foveal pit is known to be a radially asymmetric structure-it is broader in the horizontal plane compared to the vertical directions [11]. We studied 24 angular directions individually and found that percentual age and sex effects are relatively uniform across all directions. This can indicate that, despite its structural asymmetry, the fovea is a homogeneously evolving structure.
This study is not without limitations. First, using 25 B-scans under-samples the central region, which introduces a systematic overestimation and underestimation of central thicknesses and foveal slope measurements, respectively. While this is important to consider, the sensitivity analysis also revealed that the bias in most of the estimations is relatively low (<5%). We followed our standard 25-Bscan protocol as it is used in regular practice and requires a shorter scanning time, which is crucial when imaging subjects with neurodegenerative diseases. Second, we employed raster scans and interpolation to analyze the fovea radially instead of using a radial acquisition pattern. Although this may hinder an accurate reconstruction of the TRT profile in the vertical direction, we selected it because radial patterns-with an irregular sampling density-could potentially introduce a bigger bias when measuring thicknesses far from the central region or when correcting fixation errors.
We did not correct for display distortion, as neither biometric nor scanner optical information was available [49]. The errors due to such distortion are minimal for thickness measurements and small fields of view like the one used in this study [48] but have an impact on slope metrics [52]. We addressed the latter by computing the slope after flattening the retina (i.e., using the TRT). Although this procedure does not measure the actual slope seen in OCT images, it is the most common approach as it helps minimize the effect of both retinal curvature and display distortion.
Axial length is also known to influence retinal measurements. To address this, we adjusted all regression models using the scan focus, a parameter exported by the scanner that is measured when focusing the image and accounts for the refractive error of each eye [32]. Given that refractive error is closely related to axial length (R 2 > 0.72) [53], we considered the scan focus parameter as a reasonable proxy for axial length. Additionally, we relied on the lateral image scale estimation performed by the Spectralis scanner to correct the ocular magnification problem. This procedure, however, might be limited when using default corneal curvature values [32]. Finally, we did not include an interaction term between sex and age as the high inter-subject variability would diminish the statistical power to detect a potentially very small effect.

Conclusions
Most retinal layers present thinning over time that is more prominent for the GCIPL. Percentual changes in TRT and GCIPL are homogeneous even when analyzed in very small sectors. Overall, males have thicker retinal layers than females, although the differences are more evident in the inner ring. The clearest effect of age on the foveal pit is a decrease in the rim height. Male and female foveae show evident differences, with females having a shallower and broader pit. Sex and age effects are present for all angular directions. Advanced analysis of OCT images such as highly detailed thickness sectorization or radial geometrical analysis of the foveal pit can be used to enhance the description of the macula.